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Special bases of orthogonal polynomials are denned, that are suited to expansions of density and 
potential perturbations under strict particle number conservation. Particle-hole expansions of the 
density response to an arbitrary perturbation by an external field can be inverted to generate a 
mapping between density and potential. Information is obtained for derivatives of the Hohenberg- 
Kohn functional in density space. A truncation of such an information in subspaces spanned by a 
Ch , few modes is possible. Numerical examples illustrate these algorithms. 
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00 ! i. INTRODUCTION 

(N : 

J> ■ The well-known Hohenberg-Kohn (HK) density functional [1] and its finite temperature generalization by Mermin 
l/"") , [2] suffer from the absence of constructive algorithms after their respective existence theorems. The Thomas-Fermi 
(TF) approach, however, and related developments such as [3], [4], have gone a long way into creating functionals with 
^? , practical physical values. For reviews on the effectiveness of the detailed forms of the functional found empirically, 

see for instance [5] and [6] . For applications of Skyrme forces to nuclear densities, see for instance [7] and [8] . 
I/" - ) , Standard perturbation theories (particle hole hierarchy of excitations, configuration mixing, generator coordinates, 
etc.), extrapolating from well understood mean field theories, give a constructive approach to the intricacies of a true 
ground state (GS), at the well known heavy cost of calculations with many degrees of freedom. But such theories 
proved to be practical, because suitable truncations were found that restricted calculations to few modes, collective 
or not, subspaces with fewer degrees of freedom. The purpose of the present note is to attempt answering a similar 
question in the space of densities rather than the space of wave functions: are there possible truncations, is there a 
possibility to restrict the functional to a set of few density modes? 

For this we visit again the fundamentals of the HK functional F[p] in a systematic approach, based upon the 
following chain of arguments, 

i) given the full Hamiltonian, H = T + V + U, with a fixed kinetic operator T = t%, a fixed two-body potential 
operator V = X^>j w u> anc ^ a variable one-body potential operator U = assume a non degenerate, square 
normalized GS with its corresponding eigenvalue E, density p and functional F[p] = (^f\(T + V)\^f) = E— (^f\U\^); 
find the functional derivatives Sp/Su and, considering first F as a functional of u rather than p, find SF/Su, 

ii) expand such functional derivatives into suitable bases, to describe them by convenient matrices and vectors, 
hi) then invert Sp/Su to know Su/Sp, 

iv) furthermore obtain SF/Sp by eliminating Su between SF/Su and Su/Sp; further information about F might be 
obtained by integrating SF/Sp, or by comparing with phenomenological approaches, such as gradient expansions, 

v) at each stage, try a compression of the information, by a truncation of the theory to a few "density modes" . 
A preliminary question is in order, however: can this formal program be carried if particle number is conserved 

in the mean only, as occurs with Lagrange multiplier techniques? According to [5], the chemical potential, as a 
function of a continuous particle number, shows derivative discontinuities. We thus find it safer, in this paper, to 
stick to "slices" of the functional, those for fixed, integer particle numbers. We also restrict our considerations to pure 
eigenstates of H, at zero temperature. 

In Section II, we carry the program in the absence of V; the trivially soluble situation of independent fcrmions 
allows us to easily describe the mapping, p «-> u, for both infinitesimal and finite variations of u. In Section III we 
reinstate V, but use the Hartree-Fock (HF) approximation to still obtain p without excessive technical complications. 
Section IV is dedicated to a better understanding of the "tangent" mapping Sp/Su [9] when full correlations are 
present. Then Section V introduces, via a new family of orthogonal polynomials, candidates for density and potential 
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space modes, that might allow a compacted description of the functional and its p <-> u mapping. An investigation 
of the relevance of such modes is provided numerically, for a toy model of independent fermions. The numerical 
investigation is continued in Section VI, by means of a second toy model, with now correlated fermions. Section VII 
contains a discussion and a conclusion. 



II. PARTICLE-HOLE EXPANSIONS FOR INDEPENDENT FERMIONS 



Let Z be the particle number for a finite Fermion system. For simplicity we ignore discrete labels such as spins and 
isospins. Set v = 0, temporarily, the case of two-body forces being discussed later. The one-body potential function u 
is taken here as a local potential u(r). The Hamiltonian then boils down to H = Y^f=i(^i + an d its GS, assumed to 
be non degenerate, is a Slater determinant. The HK functional, in this special case, reduces to the kinetic expectation 
value, (W|T|W). 

Consider a perturbation Su of the local potential u. We find it practical to expand it in an orthonormal basis of 
functions w a (r), namely Su (r) = Yla=i w a{ r ) $ u a- According to the HK theorem [1], this basis must be orthogonal 
to a "flat potential component" wo(r) = 1. This is satisfied if we find a basis such that, Va > 0, J drw a (r) = 0. It 
is then understood that the index a will run from 1 to oo. For obvious practical reasons, however, the expansion will 
sooner or later be truncated at some finite rank of the basis. 

The perturbation Su induces a perturbation Sp of the density of the GS of H, and we find it convenient to 
expand Sp in the same orthonormal basis {wp{r)~\, namely Sp(r) = Y^pLi w p{ r ) &P@- The fact that this basis satisfies 
the constraint, V/3 > 0, J drwp(r) — 0, is very useful because Sp does not change the particle number, namely 
Sp automatically satisfies the condition J drSp(r) = 0. It is very convenient that the constraint of particle number 
conservation and that of "non constant potential variation" allow the same choice for our forthcoming basis {w a }. 

A trivial particle-hole argument then provides that perturbation <$*!> induced by Su. Let the "hole index" i = 
1, Z and "particle index" / (running from Z + 1 to oo) denote occupied and empty orbitals, respectively, with the 
corresponding single particle energies rji and r\i and orthonormal wave functions ipi and ipi. For the sake of simplicity 
in the following, such orbital wave functions tp(r) are assumed to be real in the coordinate representation. Each filled 
orbital picks a variation Stpi(r) = J2i ^i{ r ) (I\<* u \i) / {Vi ~ Vi)- Hence, 

Sp(r) = 2Y / Mr)Mr) { -^^- (1) 

This reads, when Su and Sp are expanded, 

Sp = 2 V" Vpu ^ Wa ^ § Ua (2) 

V Vi-m 

tla 

where V denotes the projection of a particle- hole product of orbitals upon the basis {w a }, 

V 0iI = J drwp{r) ^ (r) i/>/ (r) . (3) 

Note, incidentally, that particle-hole orthonormality ensures that, Mil, J dr tpi(r)tpi(r) — 0. Hence, functions w a 
expanded in the basis of particle- hole products ipiipi, a basis to be orthonormalized, automatically fulfill the requested 
condition for Su and Sp. Furthermore, positivity of the density is guaranteed as variations Sp in Eq. (1) are based on 
variations of the wave function, in particular by particle-hole admixtures to the GS determinant 
Define the matrix, 

JV 0a = 2j2v 0U ^^. (4) 

Notice also that the perturbation matrix element, (I\w a \i), coming from Su, is nothing but an integral of a three 
term product, amounting to T> a u. Notice finally that the energy denominators correspond to a propagator G = 
Q(E — QHQ)~ 1 Q, if Q is the particle-hole space projector. This operator is diagonal in the particle-hole space, 
obviously. Then, in a condensed notation, M — 2VGV, where the tilde denotes transposition between index a and 
pair index il. Consider now Sp and Su as just vectors with components Spp and Su a , respectively, sooner or later 
truncated. Since Eq. (2) reads Sp = ATSu, the HK theorem states that, under the usual condition of non degeneracy 
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for the GS of H, an inversion is possible. Namely for any Sp which leaves p + Sp in the manifold of actual densities, 
there exists a unique u + Su, provided Su does not add a constant component to u. Under such precautions, the infinite 
matrix J\f can be inverted, and the same can be expected under "reasonable" truncations of J\f. Accordingly, while 
Eq. (2) provides the functional derivative Sp/Su, one obtains the functional derivative Su/Sp, 

8u a = Y,( N ~ l ) a p 5 PfJ- ( 5 ) 



Now, that variation SE of the GS energy induced by Su is trivial. It just reads, 

5E=(V\5U\V>), (6) 
because of the stationarity of the GS energy with respect to Sty. Accordingly, for the functional under study, 

SF = SE- 5{*\U\V) = -(<5*|£7|#) - {*\U\6V) = -2<tf|CJ|Jtf) = -2{*\U G 8U\V), (7) 

hence, 

SF - -2^<i|«|/) - m)- 1 (I\6u\i), (8) 

and finally, with proper expansions, 

SF = -2 £<i|«|/) - r?,)" 1 ]T P m/ ]T (AT 1 )^ Spp . (9) 

»J a 

It may be convenient here to set from u a column vector U with components 

U u = J dru(r)Mr)ipi(r), (10) 

hence the functional derivative SF/5p reads, in a condensed, matrix and vector notation, 

5F = -UGV (vGVj 1 Sp. (11) 

This simplifies if one observes that, because of orthogonality between particle and hole orbitals, the product ipi(r) ipi{ r ) 
can be expanded in the w-basis as, 

1>i(r)Mr)='52ViipMr)- ( 12 ) 



Accordingly, 

Uu = I dr u{r) ^ V iI0 w (r) =^u l3 V iI p , (13) 

J 

with up the components of u in our special basis. Combining Eqs. (11) and (13) results in 

8F = -Y J Up(vGT>) (vGVy 1 5 Pa = - J2u a S Pa . (14) 



This avoids the transition between different bases through the matrix T>. One thus recovers the trivial result, SF — 
— J dru{r) Sp(r), but it must be kept in mind that, here, u has become a functional of p. Whether this simplification 
is made or not, this makes a set of numerical, non linear, coupled, partial differential equations relating F and p. 
The non linearity comes in particular from the orbitals and single particle energies which occur in the definition of 
N. We stress again that the vector, UV^ 1 , just makes an "a" representation of u, converted from its particle-hole 
representation U. 

A comment about Legendre transforms is here in order [9]. According to the Hellmann-Feynman theorem, SE/Su = 
p. But then, F = E — J dr u(r) p(r) is nothing but the Legendre transform of E and the primary degree of freedom is 
not u any more, but p. Note that the reasoning remains if V is reinstated. In all cases, u is recovered from SF/5p = —u. 

For Eq. (11), and its generalization if two-body forces are present, to become a tool to obtain information about 
F, dynamical models are obviously necessary. These are the subjects of several of the forthcoming sections. 
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III. TWO-BODY FORCES AND HARTREE-FOCK MODEL 



In this section, we stay with Z fermions, but reinstate in the Hamiltonian the two-body interaction Vij with the 
operator V = X)i>j=i v ij- The HK functional is (*&\(T + V)\^}. While a Slater determinant ^ was available as the 
true GS of a simpler H = T + U in the previous section, we cannot usually obtain the true GS with two-body forces 
present in a full H = T + V + U. Thus, in this section, we tolerate for ^ the HF ground state of H, with energy E , 
and furthermore assume that this HF approximation does not create degeneracies between distinct ^F's. Under this 
precaution of uniqueness, there exists an extension of the HK theorem. Indeed, in the space of Slater determinants, let 
* and be the HF GSs of H = T + V + U and H' = T+V+U', respectively, and let p{r) and p'{r) be their respective 
densities. The two Hamiltonians differ by their (local) one-body operators U = J2i=i u ( r i) an d U' — J2i=i u '( r i) 
only. Their HF GS energies Eq — (\I r |if|\E r ) and Eq = (^'\H'\^>'), non degenerate, may be equal or distinct. Now, if 
p and p 1 were equal, then the usual HK arguments, namely E' — Eq < (fy\(H' — H)\ty) = f dr [it'(r) — u(r)] p{r) and 
E Q - E' Q < (&'\(H - if') |*') = Jdr [u(r) - u'(r)] p(r), necessarily lead to p ^ p', by contradiction. 

It can be stressed here that, again because of the stationarity of the energy with respect to variations of we can 
still take advantage of Eq. (6) for the variation of the energy induced by a variation Su. This reads, with notations 
already used in the previous section, 

SE = (tf |«7|tf) . (15) 

The same holds every time we approximate the GS by means of the Rayleigh-Ritz principle in a restricted space 
of wave functions. As a general consequence, we obtain again Eq. (7), namely, SF — — 2(\^ , |f7|(5\f , ), for every such 
variational approximation of 

That variation of S^ induced by Su is slightly more complicated, in the HF case, than in the trivial case of Section 
II where H = T + U. Indeed, each filled orbital is driven by the perturbed HF equation, 

h 2 

[ Vl + S m - u{r) - Su(r)} [^(r) + <fy>;(r)] + — A r [^(r) + <^(r)] = 
z 

/ dr' v(r - r') [^(r') + <%(r')] 2 Mr) + <tyi(r)] - 

Z 

/ dr' v(r - r') [^(r 1 ) + <%(/)] [^(/) + S^(r')} [^(r) + 5^(r)\ . (16) 

3 = 1 J 

The non locality of the HF mean field, because of antisymmetrization, is written in an explicit way in the right-hand 
side above, in the coordinate representation. An equivalent form of this perturbed HF Eq. (16) is obtained if we 
retain its first order terms only and consider the particle-hole infinitesimal components Sen, again assumed here to 
be real numbers, 

(m - m) Sen - (i\Su\i) = ]T [(u\v\ij) + (ij\v\ij)] Scjj. (17) 

3 J 

Thus tpi(r) becomes ?pi(r) + ipi(r) Sen . Notice that Sr/i drops out from the calculation, as should be expected. 
Then define in particle-hole space the symmetric matrix, with antisymmetrized matrix elements of v, 

A (H)UJ) = (Vi - Vi)Si]Su - (I J\v\ij) - (Ij\v\iJ). (18) 

Here S is a Kronecker symbol and we must use pairwise indices (il) when defining the inverse A~ x to be used; this 
A~ x generalizes the propagator used in Eq. (4), and thus, 

Scu^Y,( A ^)(u)(3J) < J I<Hj)- (19) 

3 J 

This leads to the variation Sp, and the analog of Eq. (2) reads, 



5 Pfi = 2 J2 V fjU (A- 1 )^^ (J\w a \j)5u a , (20) 

(U)(jJ)a 
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where the overlap matrix V is the same as defined in the previous section. Actually, this boils down to the even 
simpler formula, in matrix and vector notations, 

Sp = 2VA- 1 V5u, (21) 

where the tilde again denotes transposition of that connection T> between the particle-hole products , 4'i{r)ipi{r) and 
their rearrangement into an orthonormal basis {w a (r)}. Note, incidentally, that, if v = 0, the matrix T = 2D A^ 1 D 
boils down to the matrix A/", which is obviously negative semidcfinite. We even expect Af to be negative definite. The 
same is expected for D A^ 1 T>. The stability of our HF solutions is assumed as long, at least, as v is a weak enough 
interaction, and this "definiteness" of T is intuitively most likely. 

In the following, we shall also need the inverse of T . The final result for the variation of F reads, 



5F=-UA- 1 f> VA^V 



-l 



Sp, (22) 



and, like in Section II, this expression, in a transparent notation, reads 8F = — ~^2 a u a 5p a . For obvious reasons 
of numerical convergence, the number of needed w a states must be large enough to overlap a sufficient number of 
particle-hole components of 5p. But as will be found in the coming numerical applications, a surprisingly small number 
of w a states might sometimes be sufficient. 

IV. SYMMETRY OF THE DENSITY-POTENTIAL MAPPING IN GENERAL 

With the exact ground energy E and exact GS W of a full H = T + V + U, and Q = 1 - |*)(*| the projector out 
of the Brillouin-Wigncr perturbation theory gives the exact result for first order functional derivatives, 

n— 1 

We assume here that a resolution of the identity with real numbers and reasonable truncations, convenient for numerics, 
are available. The sum over excited states \I/„ includes integrals over the continuum, if necessary. Let us single out 
the first of our identical particles and integrate out all the other ones, to define the following transition densities, 

0„(r) = J dr 2 dr 3 ... dr z ^ n (r,r 2 ,r 3 , ...,r z )^(r,r 2 ,r 3 , ...,r z ). (24) 

Notice that, from its very definition, 0„ integrates out to 0, namely J dr n (r) = 0. Hence n can be represented in 
the w-basis without any loss of information. 



Since SU = J2i=i <^ u ( r i) is a symmetric operator, it is clear that Eq. (23) also reads, 

m = g ^ ZjdrSu(r)e n (r) = ^ ^ Z J ^(r) 6 m (r) ^ 



E — E n ^ E — E n 

n=l n a 



(25) 



where we have again expanded 5u in the basis {w}. There pops out a matrix, 

D an = Z J drw a {r)Q n (r), (26) 

as a generalization of the matrix T> a u. 

Now, by definition, the density of the GS is, 

p{r) = Z j dr 2 dr 3 ...dr z [*(r, r 2 , r 3 , r z )f , (27) 

and its variation is, 

Sp(r) = 2Z j dr 2 dr 3 ...dr z *(r, r 2 ,r 3 , ...,r z ) Sf(r, r 2 ,r 3 , ...,r z ). (28) 
This becomes, if one replaces (5^ by its expression, Eq. (25), 
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E — E n 



An expansion of dp in the {w} basis gives its coordinates, 



5 P0 = 2'£D i 



1 



pn 



E — E n 



D a n 5u a7 



(29) 



(30) 



hence, in an obvious notation, a symmetric "flexibility" matrix T = 2DG D connecting Su and 5 p. In hindsight, the 
symmetry of T (and of its approximations under the Raylcigh-Ritz variational principle) is straightforward. Indeed, 
since u = — 5F/Sp, then Su a /Spp = — S 2 F/(5p a Spp). All denominators E — E n being negative definite, the negative 
definite nature of this exact T is also transparent. 

We conclude this section on the general case with explicit expressions for 5F/5p and S 2 F/(5pSp'). With Eqs. (25), 
(26) the general form for 5F, see Eq. (7), reads 

SF= -2{Sf\U\8*) = -2Y J {^\U\* n ){E-E n )- 1 D an 5u a = -2^JA n (E - E n )~ 1 D an Su a . (31) 

na na 

Here the numbers 



U n = J dr<*| 5>(r0 |*„) = Z J 



dru(r) 9„(r), 



(32) 



are now the components of u in the space of transition densities, generalizing Eq. (10) for states containing 
correlations. Upon inverting Eq. (30) we find, as a generalization of Eq. (22), 



SF = -UGD 



DGD 



Sp. 



This simplifies if we expand 



Then 



® n (r) = Z- 1 Y J D n fiW f} (r). 





U n = dru(r) D n/3 wp(r) = ^2 up D n p, 
J fi fi 

with again the components of u in our special w-basis, up = J dr u(r) wp(r). Accordingly, 



UGD 



1 



= y^jUfi Dfin Gnn D nl — - ^ Up Tp 1 , 



fin 



hence finally, as expected, 



SF^-^upTp^ [T 1 ] rya 5p a = -^2u a 5p a , 



(33) 



(34) 



(35) 



(36) 



(37) 



as before in sections II and III. Similarly, the second derivative of F is found directly from the inverse of Eq. (30), 



5*F/(5 Pa 5pp) = {T- l ) af: 



(2DGI)) 



afi 



Hence, from Eqs. (33) and (38), 



SF/S Pa = 2 J2 (UGD) 5 2 F/{5 P p5 Pa ), 



(38) 



(39) 



another useful equation to test phenomenological functional F[p], by calculating quantities such as UGD from 
microscopic wave functions and energies for simple systems. 
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V. ONE DIMENSIONAL TOY MODEL, SPECIAL POLYNOMIALS 



Assume that r is just one dimensional, running from — oo to oo. Define t = —d 2 / (2dr 2 ), with a nucleon mass m = 1, 
h = 1 and p = —id I dr. In the present section, we are first interested in the Hamiltonian Hq = X)^=i(-Pi/2 + T 1l^)\ 
it is not a bad approximation to most shell model Hamiltonians, whether one considers one-body potentials only or 
HF solutions to problems with two-body potentials as well. A trivial scaling of coordinates and momenta allows us to 
reduce to the case, 10 = 1, any situation, (pf/2 + w 2 r 2 /2), where the physical spring constant would be different. 

Then we shall consider the functional F(p) = (^>\Ho\^f) for a family of additional one-body potentials u, with the 
corresponding GS density p(r) of 

z 

H = Y J \U+r 2 l /2 + u{r i )]. (40) 

i=l 

Set temporarily u — 0, namely consider H and its GS density pa(r). Since u) = 1, which will be understood from 
now on, both initial particle and hole orbitals ipk{r) are just trivial products ipk of a Hermite polynomial, a common 
Gaussian and a suitable normalization, (pk( r ) = ?r~3 e~2 r Pk(r). For the sake of illustration, we list here the first 
five Hermite polynomials, with their coefficients adjusted for orthonormalization, 

P = 1, Pi = 25 r , P 2 = 2"5 (2r 2 - 1), P 3 = 3"5 r (2r 2 - 3), P 4 = 2" 1 6 - ' (4r 4 - 12r 2 + 3). (41) 

Particle-hole products, <fi{r) ipi(r) 7 make, in turn, just polynomials again, now multiplied by e~ r . To build our 
basis, {u>Q,(r)}, it is tempting to orthonormalize the set {ipiipi} containing that Gaussian, e~ r , and recover forms 
23 tp k (r^/2), with Hermite polynomials again, compressed by the obvious \/2 for their argument r, because of the new 
factor, e~ r . This is correct for odd parity functions. But, for even parity ones, the constraint, / drw a (r) = 0, would 
be violated. Hence, out of each even function, 2« ^kirVty, k > 0, we subtract a term proportional to 2^ipo(r\/2), 
letting the subtraction cancel the integral, . (Alternately, we considered all elementary functions r 2k e~ r .) Then 
we use a Gram-Schmidt process to reorthonormalize such subtracted states. Notice that the subtraction cancels out 
the polynomial state of degree zero, and therefore the transformation from Hermite polynomials to this new set of 
orthonormal polynomials is not unitary, but only isometric, with "defect index" 1. In other words, our basis has 
codimension 1. This is also clear from the degree 2 of the lowest member of the new even basis. For an illustration, 
we list the first four even states obtained, 

w 2 (r) = 2 23 (2r 2 - 1)/Vs tt^ e'^ , 
W4 (r) = 23 (8r 4 - 14r 2 + l)/VTE tt^ e~ r * , 
w 6 (r) = (32 r 6 - 128 r 4 + 94r 2 - ll) / (2iVl0^j tt^ e~ r \ 

w s (r) = (I28r 8 - 928r 6 + 1752r 4 - 906r 2 + 39) / (9 23^35^ W e -'' 2 . (42) 

For the sake of comparison with Hermite polynomials, which rather go with a factor e~^ r , we perform the trans- 
formation, r — ► rj\pl on Eqs. (42) and multiply the results by a factor 2 _ 3 to retain their (ortho)normalization. 
Discarding norm coefficients from the resulting polynomials we get, 

Q 2 = r 2 -l, 
Q 4 = 2r 4 -7r 2 + l, 
Q e = 4r 6 -32r 4 + 47r 2 - 11, 
Q 8 = 8 r 8 - 116 r 6 + 438 r 4 - 453 r 2 + 39. (43) 
2 

But, as already noticed, products ipi ipj carry a factor e~ r and we find it natural, in the following, to stick to those 
polynomials trivially derived from Eqs. (42), 

r 2 (r) =2(2r 2 -l)/V3, 
r 4 (r) = 25 (8r 4 - 14r 2 + l)/\/l5, 
T 6 (r) = (32 r 6 - 128r 4 + 94r 2 - ll) /V210, 

r 8 (r) = (l28r 8 - 928r 6 + 1752r 4 - 906r 2 + 39)/(l8\/35) , (44) 
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and so on. We generated such polynomials up to decree 100 and will send them to interested readers. Such polyno- 
mials are orthonormal under the metric weight, e~ 2r They must be completed by odd Hermite polynomials, 
^2fc+i(f\/2), suitably adjusted for the same metric. Hence, for instance, 



(45a) 
(45b) 
(45c) 
(45d) 

More technicalities on such polynomials T and related polynomials can be found in [10]. 

It is then trivial to calculate both even and odd blocks, respectively, of the initial matrix T>, see Eq. (3), according 
to the parity of the subscript of w and its associated polynomial T. With due normalizations, this reads, 



ri(r) = 2r, 
r 3 (r) =r(4r 2 -3) y/2/3 , 
T 5 (r) =r(16r 4 -40r 2 + 15)/\/30, 
T 7 (r) = r (64r 6 - 336r 4 + 420r 2 - 105)/(6 V%5) . 



Van = 



dr 



T (r)e- r (2/tt)* P^r) e~^ r Tr"* P/(r) e _ * r 



(46) 



For instance, if the hole label is restricted to i = 0, and the particle label I runs from 1 to 3, the sets of non vanishing 
odd, respectively even, matrix elements boil down to, 



V 



101 



2~* n~*, V 



103 



V3/ U(2tt) 1 



/4 



V. 



303 



2 4 7T 4 , 



V 



202 



V3 T 



4 7T 4 . 



(47) 



We found it useful to precalculate and store such initial matrix elements V for the particle index I running up to 
100 and the a index running up to 100 also. This fastens generic calculations of T> when u becomes finite, as one 
represents (t + r 2 /2 + u) by a matrix on the oscillator basis, diagonalizes it with eigenvalues r)k and orthonormal 
eigenvectors Xik, and finally expands orbitals of both holes and particles as tpkir) = ^2 e Xgk ¥v( r ) m the same basis. 

In [10] we set Z = 4, considered u to be an infinitesimal 5u in the neighborhood of u = 0, then calculated and 
diagonalized the functional derivative J\f — 8p/Su. The eigenvectors of Af defined density and potential infinitesimal 
perturbations having the same shapes. Now we set again Z = 4, but are rather interested in cases where u is finite. 
We are concerned in particular with the mapping between u and p, in that representation provided by the "modes" 
w a . Truncations at a maximum degree N are necessary. The finite expansion, 



N 



i(r) = ^2 u aW a (r), 



(48) 



a=l 



defines those processed perturbations u. Given u, it is trivial to diagonalize H with a good numerical accuracy and 
obtain p. Then it is easy to obtain "coordinates in density space," 



p a = 



drw a (r) [p(r) - p (r) 



(49) 



The harmonic potential, r 2 /2, serves here as the origin in potential space, and the origin in density space is the 
corresponding density po- We show in Figures 1 and 2, respectively, a grid of values {^2,^4} in potential space and 
its image grid of density coordinates {p 2 , pi}. Dots at grid corners help matching the object and the image. 
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FIG. 1. Grid of parameters ui, 114 for the potential u — ui 102 + U4 W4 used in the toy model. 
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FIG. 2. Density space image, projected onto the p2,pi plane, of the grid of potentials of Fig. 1. 
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FIG. 3. Same as Fig. 2, but now the image grid is projected onto the pQ, Ps plane. 




FIG. 4. Toy model HK functional in a {p2,pi} frame. Note small deviations from paraboloid. 

In this calculation, all coordinates u a have been set to vanish, except u-i and M4, but it must be stressed that 
the resulting density variation, p — po, has non vanishing coordinates pe,pg,,... besides pi and p±. Such additional 
coordinates are small, but not very small, as shown by the grid for pe,Ps in Figure 3. Qualitatively, if u contains one 
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mode w a only, then pp tends to decrease when \a — (3\ increases. But this is likely to be valid for small enough ti's 
only, in a linear response regime. Curvature effects, evidenced by Figs. 2 and 3, must be expected further. 

Of interest are plots of F in p-space. If u has two components u 2 ,u 4 only, assume that F is a function of p 2 ,P4 
only. Then a gradient VF(p) can be observed directly. In Figure 4, the 3D plot of F shows slight deviations from a 
traditional paraboloid. This is even more visible in Figure 5, showing the vector field {u 2l ua\{p 2i Pa), namely — VF. 
The field, read from Figs. 2 and 1, focuses towards the origin in p-space, but with clear distortions. We know that 
the field has a vanishing curl; it can be integrated back into F. 

It is also trivial to create an approximate F in the following way: i) assume indeed that F depends only on p 2 and 
Pi for this toy model, ii) take a few exact (numerical, actually) values of F at random points taken from the partner 
grids shown in Figs. 1 and 2, iii) set a simple parametric ansatz such as, 

F apP F oa + F w p 2 + F01 p 4 + {F 2a pl + 2 F n p 2 p 4 + F 02 pfj /2 + (F 30 pi + 3 F 21 p 2 2 pi + i F 12 p 2 p\ + F 03 pi) /6 + 
(F 40 P \ + 4 F31 p^p 4 + 6 F 22 p\ pi + F 13 p 2 p\ + F 04 pf) /24 , (50) 

and, finally, iv) least square fit the "exact" values selected at step ii). There are here 15 parameters and it is reasonable 
to select typically about twice as many exact values for the least square fit. The following result, 

F app ~ 8.0005176 - .0026263 p 2 + 3.7994711 p\ - .9776987 p\ - .0208120 p\ + .0034190 p 4 + 
.9487531 p 2 p 4 - .2536905 p\ p 4 - .5398019 p\ p 4 + 3.7854603 p\ + .9596720 p 2 p\ + 

.0368183 p\p\ + .0172911 pi - 3.0193024 p 2 pi - .2356997 p\, (51) 

comes from fitting 26 values for Z = 4. Figure 6 shows several resulting contours, the smallest of which locates the 
minimum of F app very slightly only away from the origin, that is the true minimum by the very construction of the toy 
model. The value of the functional at the minimum turns out to be 8.0005165, instead of strictly 8. This toy numerical 
exercise demonstrates the possibility of contracting the description of the functional to few degrees of freedom. 



rho4 




rho2 



FIG. 5. Same as Fig. 4. The lines represent — VF at points {p2,P4} shown by dots. Note deviations from radial pattern. 
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FIG. 6. Contours for F c 



VI. CORRELATIONS, FROM ANOTHER TOY MODEL 

In the previous sections we skirted around the difficulty of obtaining a GS with true correlations. Now we shall 
mix several Slater determinants, each made of Z = 4 harmonic oscillator, one dimensional orbitals taken from 
ho = \ {p 2 + r 2 ). The Hamiltonian is a complete one, 

4 

-V a J2 [Sin-rj-RJ+Sin-rj+Ra)]. (52) 

i>j=l 

This contact, finite range attraction between particles is expected to create a reasonable amount of correlations and 
was chosen to allow an easy precalculation and tabulation of matrix elements of v. Such matrix elements are again 
understood to be antisymmetrized. 

The first Slater determinant, $o> i n the mixture contains the lowest Z orbitals of the harmonic oscillator. It is 
expected to make the dominant component of the configuration mixture \&, as we shall keep 112, U4 within the grid seen 
in Fig. 1, and also the strength V a moderate. Let £j and ipj, i,j = 1, Z be the orbitals of two Slater determinants S 
and $, respectively. Define the cofactors CV, and double cofactors Cikji of the determinant of scalar products (CilVj)- 
Such cofactors are very simple in the present orthogonal basis of orbitals, obviously. Then the matrix elements needed 
for the Hamiltonian matrix and the calculation of the HK functional read, 

(s|(jr + = J2 c ^ telC* + «)bfc>. = \ E &£M<pm) • ( 53 ) 

ik ijkl 

With V a = 5 and R a = 1 the grids shown by Figures 7 and 8 come from a calculation with a single particle basis 
made of the first 10 harmonic oscillator orbitals and a corresponding Slater basis of 61 states made of $0 and all 
positive parity one-particle-one-hole and two-particle- two- hole determinants built upon $o- Because of V the centers 
of the density grids are not at the origin defined by $0, obviously. Indeed, for U2 = M4 = 0, this calculation gives 
{p2, Pi, ...P14} = {—-14, —.36, —.43, —.06, .04, —.01, .005} as the coordinates of the ground state density shift p — pa 
and the ground state energy is E = —6.7, significantly down from (<&o|-£fo|3>o) = 8. 



+ u 2 w 2 (ri) + U4 W4,(n) 
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FIG. 7. Second toy model: influence of the two-body force V on the p2,p4 image of the grid of Fig. 1. Compare with Fig. 2. 




FIG. 8. Same as Fig. 7, but projection of the image grid into the P6,ps plane. Compare with Fig. 3. 
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FIG. 9. HK functional of the second toy model in {p2,Pi\ frame. 



Had we taken as the origin in density space the density of the HF solution for u = 0, different drifts of grid centers 
would have been observed. Such new drifts are likely to make better signals of true correlations in ^. For the sake of 
comparison between the first and the second toy models we kept po as a reference, but we have a direct access to the 
amount of true correlations: it is easy here to calculate the density matrix p, the diagonal of which gives p. Here, for 
u = 0, the trace of p — p 2 is of order 5%, a reasonable amount. At that grid corner, u-2 — U4 — —2.2, the trace even 
reaches 8%. It can be concluded that this second toy model does create correlations. 

While P2 and P4 both vary by ~ .3 across their grid and pw,Pi2--- can be neglected, the variation of pe and pa 
across the grid is of order ~ .1, which is not so small. For the sake of comparison with Fig. 4 we now show in Figure 
9 a plot of the HK functional in terms of {p2, Pa} again, but it will not be not forgotten in our agenda to find those 
two orthonormal combinations {02! 04} of P2, pi, Pe and pg, which allow the "flattest" projection of the grid. It is 
clear that, in the spirit of Eqs. (50) and (51), the best parametrization of F should now be in terms of {02, 04}- In 
any case, from Fig. 9, the minimum of F occurs at {p2,Pi} = { — -14, —.36}, the grid center, as should be. Fig. 9 also 
shows deviations from a paraboloid clearly stronger than those of Fig. 4. 

VII. SUMMARY, DISCUSSION AND CONCLUSION 

Particle number conservation is essential to that key element in the proof of the HK theorem, the one-to-one 
correspondence between density and external potential. However, variations of the potential should not be trivial: 
they should differ from constants. We treated both constraints of i) matter conservation and ii) non triviality of 
potential variations on the same footing: a vanishing average for both the potential and the density variations. 

This constraint of vanishing average was implemented by means of a new family of orthogonal polynomials; hence 
appeared a set of modes in both the density and the potential spaces. We proved, numerically with toy models, 
that such modes might have a physical meaning, on two counts, i) converse linear responses Sp/Su, du/Sp might 
be reasonably simple when described in terms of such modes, and ii) the HK functional itself might be practically 
truncated into projections into subspaces spanned by a few modes. 

We have not discussed in this paper the constraints of positivity of the density, but it is clear that, within an 
algebra of polynomials such as ours, positivity conditions are not too difficult to implement. We have not discussed 
either more subtle constraints related to the Sobolev nature of the topological spaces available for densities. For this 
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question, we refer to [11], [12], [13]. It can be stressed again that the compatibility of truncations of densities into a 
finite number of "polynomial modes" with such fine constraints can easily be tested. 

An ultimate goal would be to create a constructive theory of the HK functional. The functional differential equation, 
5F/5p = —u[p], cannot be integrated in the density space as long as u[p] is not known accurately enough. Because of 
our detour through many-body perturbation theories we are clearly far from the goal, but this work gives a frame in 
which the task should become easier. The detour might allow a compression of the needed information through, for 
instance, the parametrization of a limited set of matrix elements of our matrices T> for mean field approximations or D 
in the correlated cases. Our main results are i) the existence of those special, orthogonal constrained polynomials r a 
and associated modes w a which design convenient sets of coordinates and convenient parametrizations of F, SF/Sp, 
S 2 F/ (SpSp 1 ), ... etc., ii) the explicit relation we showed between these modes and the traditional perturbation theories 
used in the many-body problem and iii) the likely possibility of truncated descriptions and accurate parametrizations. 

Actually, in the context of extended systems, the idea of density waves as important modes of the system has always 
been present. There remains to be seen, obviously, if our modes can be generalized to infinite systems. It also remains 
to be seen whether, for finite or infinite systems, our truncations are always justified, whether infinite resummations 
are possible, whether collective degrees of freedom are present, or absent, because of our new representation. Also, 
because of our choice of a Gaussian weight for the new family {T}, the present results arc better meaningful if restricted 
to nuclei. A generalization to atoms and molecules obviously demands other weights for the constrained polynomials. 

The usual perturbation theories have their hierarchy of modes in many-body space, most often a hierarchy of particle- 
hole components. Our approach, tuned to the one-body nature of the density functional, replaces the particle-holes 
by other modes, in a transparent way. 

It is a pleasure to thank Y. Abe, R. Balian and B. Eynard for stimulating discussions. 
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